#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "common.h"
#include <list>
#include <omp.h>

using namespace std;
double cutoff__ = cutoff;

//
//  benchmarking program
//
int main( int argc, char **argv )
{    

    //omp_set_num_threads(24);

    if( find_option( argc, argv, "-h" ) >= 0 )
    {
        printf( "Options:\n" );
        printf( "-h to see this help\n" );
        printf( "-n <int> to set the number of particles\n" );
        printf( "-o <filename> to specify the output file name\n" );
        return 0;
    }

    int n = read_int( argc, argv, "-n", 1000 );

    char *savename = read_string( argc, argv, "-o", NULL );
    
    FILE *fsave = savename ? fopen( savename, "w" ) : NULL;
    particle_t *particles = (particle_t*) malloc( n * sizeof(particle_t) );
    set_size( n );
    init_particles( n, particles );

    // Initialize Boxes! woo!
    double size = 0.0;
    size = sqrt(density*n);
    int numBoxes =(int) ceil(size/cutoff__);
    list<int> boxes[numBoxes][numBoxes]; // [x][y] 
    omp_lock_t box_locks[numBoxes][numBoxes];  //[x][y]

    double simulation_time = read_timer( );

    #pragma omp parallel
    {
        #pragma omp for 
        for (int i=0; i<numBoxes*numBoxes; i++) {
                omp_init_lock(&(box_locks[i/numBoxes][i%numBoxes]));
        }
    
        //  simulate a number of time steps
        for( int step = 0; step < NSTEPS; step++ )
        {
          
            // update boxes 
            #pragma omp for 
            for( int i = 0; i < n; i++){
                int x = (int)(particles[i].x/cutoff__);
                int y = (int)(particles[i].y/cutoff__);
                omp_set_lock( &(box_locks[x][y]) );
    
                // critical section for our stuff.
                boxes[x][y].push_back(i);
    
                omp_unset_lock( &(box_locks[x][y]) );
            }
    
            //  compute forces
            #pragma omp for
            for( int i = 0; i < n; i++ )
            {
                particles[i].ax = particles[i].ay = 0;
                int currBox_x = (int) (particles[i].x/cutoff__);
                int currBox_y = (int) (particles[i].y/cutoff__);            
    
                list<int> neighbors;
                for(int I = max(0,1-currBox_x); 
                        I < 3 && currBox_x-1+I<numBoxes; 
                        I++) {
                    for(int J = max(0,1-currBox_y); 
                            J < 3 && currBox_y-1+J<numBoxes; 
                            J++){
                    list<int>::iterator iter;
                    for(iter=boxes[currBox_x-1+I][currBox_y-1+J].begin();
                            iter!=boxes[currBox_x-1+I][currBox_y-1+J].end();
                            ++iter) {
                            // inserting all the neighboring particles to the
                            // neighbors list
                            neighbors.push_back(*iter);
                        }
                    }
                }
                // we sort to achieve perfect correctness with original
                // serial implementation TA's gave us 
                neighbors.sort();
                
                for (list<int>::iterator iter=neighbors.begin(); iter != neighbors.end(); ++iter) {
                    apply_force( particles[i], particles[*iter]);
                }
     
            }
    
            
            // clean boxes
            #pragma omp for 
            for( int i = 0; i < n; i++){
    
                int x = (int)(particles[i].x/cutoff__);
                int y = (int)(particles[i].y/cutoff__);
                omp_set_lock( &(box_locks[x][y]) );
    
                if (! boxes[x][y].empty()) {
                    boxes[x][y].clear();
                }
    
                omp_unset_lock( &(box_locks[x][y]) );
            }
    
            //  move particles
            #pragma omp for
            for( int i = 0; i < n; i++ ) 
                move( particles[i] );
    
            //  save if necessary
            #pragma omp master 
            if( fsave && (step%SAVEFREQ) == 0 )
                save( fsave, n, particles );
        }
    }
    simulation_time = read_timer( ) - simulation_time;
    
    printf( "n = %d, simulation time = %g seconds\n", n, simulation_time );
    
    free( particles );
    if( fsave )
        fclose( fsave );
    
    return 0;
}
